function neuronalDensityCalcs
% neuronalDensityCalcs

% gs Oct 2004
% ------------------------------

% neuronal densities
D = 50000; % in neurons/mm3
D = D / 1000^3; % convert to neurons/um3

% vector for diameter of volume
X = 0:30; % in um

% volumes for given diameters
V = (4/3) * pi* X.^3; % spherical model; in um3
% V = X.^3; % cube volume, in um3

% number of neurons in excitation volume
N = D * V;
V_oneCell = 1/D;
% diam_oneCell = (V_oneCell)^(1/3);
diam_oneCell = ((1/((4/3) * pi)) * V_oneCell)^(1/3);

V_twoCells = 2/D;
% diam_twoCells = (V_twoCells)^(1/3);
diam_twoCells = ((1/((4/3) * pi)) * V_twoCells)^(1/3);

% mean distance between somata as a function of neuronal density
% the value is approximated as the length of the side of the volume that holds one cell, 
% which is the inverse of the density; e.g. 10 cells/m3 ==> (X m3/cell)^(1/3)
Dvec = [20000:10000:200000];
Dvec = Dvec / 1000^3;
Xvec = (Dvec.^(-1)).^(1/3);
D_Keller = [70000 110000 50000];
D_Keller = D_Keller / 1000^3;
X_Keller = (D_Keller.^(-1)).^(1/3);

% figure -----------------------------------
figure(80);
clf;
set(gcf, 'Color', 'w', 'Position', [156   580   582   528]);

LineWidth = 1;

subplot(2,2,1);
hplt = plot(X, V, 'b-');
set(hplt, 'LineWidth', LineWidth);
xlabel('Diameter of sphere (um)');
ylabel('Volume of sphere (um3)');
pbaspect([1 1 1]);

subplot(2,2,2);
hplt = plot(X, N, 'b-');
set(hplt, 'LineWidth', LineWidth);
xlabel('Diameter of sphere (um)');
ylabel('Neurons per sphere');
pbaspect([1 1 1]);
% lines corresponding to 1 cell
hline(1) = line([0 diam_oneCell], [1 1]); % horizontal line
hline(2) = line([diam_oneCell diam_oneCell], [0 1]); % vertical line
set(hline, 'Color', 'g', 'LineStyle', '-', 'Marker', 'none', 'LineWidth', LineWidth);
% lines corresponding to 2 cells
hline(1) = line([0 diam_twoCells], [2 2]); % horizontal line
hline(2) = line([diam_twoCells diam_twoCells], [0 2]); % vertical line
set(hline, 'Color', 'r', 'LineStyle', '-', 'Marker', 'none', 'LineWidth', LineWidth);

subplot(2,2,3);
hplt = plot(Dvec, Xvec, 'b-');
set(hplt, 'LineWidth', LineWidth);
xlabel('Neuron density (cells/um3)');
ylabel('~Inter-soma distance (um)');
hold on;
hplt = plot(D_Keller, X_Keller, 'ro');
pbaspect([1 1 1]);

% sediment:
% % V = X.^3;
% Y = D * X.^3;
